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Models ol relativistic heavy ion collisions typically involve both a hydrodynamic module to de- 
scribe the high density liquid-like phase and a Boltzmann module to simulate the low density break- 
up phase which is gas-like. Coupling the prescriptions is more complicated for viscous prescriptions 
if one wants to maintain continuity of the entire stress-energy tensor and currents. Derivations for 
the viscosity for a gas are reviewed, which then lead to expressions for changes in the phase space 
occupation based on simple relaxation-time pictures of viscosity. These expressions are shown to 
consistently reproduce the non-equilibrium components of the stress-energy tensor. An algorithm 
for generating a Monte Carlo sampling of particles with which to initiate the Boltzmann calculations 
is also presented. 

I. INTRODUCTION AND BASIC THEORY 

In the past few years, the community modeling heavy ion collisions has made tremendous progress incorpo- 
rating viscous effects into hydrodynamic models. However, even viscous models are based on an assumption 
of short mean free paths and become invalid near the breakup stage when densities become low enough that 
particles can traverse a significant fraction of the reaction volume. Fortunately, since these densities are low, the 
system can be treated as a hadronic gas undergoing binary collisions and Boltzmann treatments are justified. 
This paper focuses on issues related to the interface between the two treatments. For transitions from ideal 
hydrodynamics, the transition is treated by generating particles consistent with a thermal distribution as a 
source term for the Boltzmann transport. Here, we present and test a detailed prescription for consistently 
incorporating viscous effects into such an interface. We also address the question of whether bulk viscosity plays 
a role in a low density hadronic gas. 

From a phenomenological point of view, an understanding of these issues is crucial before the viscosity of the 
so-called "perfect fluid" created in heavy ion collisions can be properly ascertained. Most of the elliptic flow 
observed in the flnal state is generated in the earlier stages of the collision [1] whic can be undestood with ideal 
hydrodynamics. However, there are substantial differences between ideal and viscous hydrodynamics. Only 50% 
of the difference of elliptic flow between ideal and viscous hydrodynamics comes from the reduction of the second 
component of the flow gradient by viscous forces. The rest comes from a modiflcation of the Cooper-Frye |2] 
formula at freeze-out [3l |4j . While this modiflcation is generally demanded by energy-momentum conservation 
in the presence of a pre-existing flow gradient, non-equilibrium effects at freeze-out are now known [S] to be 
capable of modifying the V2 resulting from a given flow gradient by parametrically non-negligible factors. Thus, 
the detailed dynamics of the decoupling of the fluid into particles, and their subsequent interactions, introduce 
a 50% systematic effect in any determination of 77/s from comparison of hydrodynamic simulations with data, 
and needs to be properly accounted for. 

A hydrodynamic model provides the fluid velocity u", the rest-frame energy density e and charge densities 
p. Equivalently, e and p can be represented by a temperature T and chemical potentials p. Assuming the 
transition density is sufficiently low to warrant treatment as a thermalized gas, the phase space densities of the 
various species are straight-forward to generate. Viscous treatments introduce new quantities describing the 
spatial components of the stress-energy tensor in the fluid frame. Whereas the components are Tij = PSij at 
equilibrium, where P{e,p) is the pressure, Tij has additional terms in a viscous treatment. In Navier-Stokes 
treatments, the additional terms are uniquely determined by the local velocity gradient, 

T^, = F<5,,+7rj;)+^W%, (1) 
^If = ~V [^^VJ + d,v, - (2/3)(5„-V • v] , tt^^^ = -CV • v, 

where rj and are the shear and bulk viscosities respectively. In Israel-Stewart approaches the shear and bulk 
corrections, tt^"^ = T^j — (Tr Tij)Sij/3 and tt*^''' = (Tr Tij)/3 — P, are treated as dynamic objects which relax 
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toward Navier-Stokes values [StilOj. A consistent interface needs to maintain the continuity of tt^^^ and 7r('') as 
well as the densities and collective velocities. 

In the next two sections, we present calculations of the bulk and shear viscosities for low density hadronic 
gases in terms of collision times. Although the shear calculation is a familiar result [Hldlj, the bulk calculation 
is more complicated, but similar calculations can also be found in the literature |13j . Sophisticated calculations 
of the shear and bulk viscosity would account for differing relaxation times for different species or for particles 
with different momentum 14 , something which could give rather non-trivial results close to the Hagedorn 
temperature |15j . By using simulations, calculations of the shear viscosity can more accurately incorporate 
effects from memory being transferred to colliding partners [1^]. The calculations of the next section are not 
meant to be competitive with these more realistic treatments. Instead, these simple relaxation-time based 
calculations are included to provide background and context for the algorithms presented here to alter phase- 
space densities so that they are microscopically consistent with the relaxation time picture. By adjusting the 
relaxation time, one can reproduce any desired viscosity. 

Although the bulk viscosity is zero for a gas of massless particles, or for a gas of non-relativistic particles, a 
small and nearly negligible bulk viscosity ensues for a gas of semi-relativistic particles. For a gas with a mixture 
of non-relativistic and highly relativistic particles, e.g. pious and nucleons, one finds a larger, but still small, 
bulk viscosity. Identical results are derived from the perspective of Kubo relations, and from a more microscopic 
picture based on the evolving phase space density in the presence of velocity gradients. The next two sections 
present results for the bulk and shear viscosities respectively, and the subsequent section describes a consistent 
design of an interface between hydrodynamic and Boltzmann models. The final section summarizes the findings. 



II. THE BULK VISCOSITY OF A DILUTE GAS 



The bulk viscosity will be non-zero whenever the trace of the stress-energy tensor, Tu/S, can differ from the 
equilibrium pressure, -P(e, p). The inability to maintain equilibrium is assumed to derive from rapidly changing 
densities, i.e., V • w ^ 0. For a non-interacting gas of non-interacting particles, 

^ = E(^/ d'p E{p,m,)Mp), (2) 

where fe{p) is the phase space density for particles of species £. For massless particles, E{p,m£) = \p\, and, as 
long as the gas is isotropic (/(p) depends on just \p\ and not on the angle of p with a preferred direction), one 
finds 

(3) 

i 

regardless of the form for /(p). Thus, there is no bulk viscosity for massless gases as the effective pressure is set 
as a function of the energy density, and is independent of kinetic thermalization. Similarly for non-relativistic 
gases, E{p,mf) = m + |pp/2m, and 

Y^Tu/3=lie-J2^ePi). (4) 

Again, this depends only on the energy and particle densities, and is independent of whether fi{p) is in an 
equilibrated form and the gas has zero bulk viscosity. A gas of semi-relativistic particles can have Tn/Z ^ P 
by distributing the energy more or less amongst the higher or lower momentum modes. Such a gas with T ~ m 
can support a small, but negligible for heavy ion collisions, bulk viscosity. The results above are common 
knowledge in the relativistic heavy ion community. What is less known is that a modest bulk viscosity ensues 
for a gas with a mixture of non-relativistic and relativistic particles. This will be illustrated below. 

For interacting systems, significant bulk viscosities can arise from several non-equilibrium effects [171 118) . 
Near a phase transition, where the vacuum expectations of fields may change suddenly, the finite equilibration 
time can lead to a peak in the bulk viscosity near Tc- Hints of this behavior can be seen in lattice results [El [20] 
and AdS/CFT calculations HHI^, and the peak in viscosity could well have a significant effect on subsequent 
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dynamics |23H25| . Long-range correlations from finite-range interactions or structures, e.g. polymers, experience 
frictional heating in an isotropic expansion and thus have bulk viscosities. Non-equilibrium chemistry can also 
lead to a bulk viscosity, but only if one compares X^i'^ii/^ to the pressure calculated with the equilibrium 
density. If the densities of various species are calculated dynamically, and if non-equilibrium concentrations are 
then used to calculate the pressure, there is no need to add a bulk viscosity to account for the non-equilibrium 
pressure. Similarly, if non-equilibrium mean fields are modeled dynamically |26j the bulk viscosity associated 
with non-equilibrium fields can be ignored. In fact, incorporating the effects of non-equilibrium chemistry or 
fields through a bulk viscosity is potentially clumsy as the effects can easily become non-linear. For this study, 
we ignore the effects mentioned in this paragraph and consider only the case of a dilute gas where the viscosities 
are solely due to kinetic non-equilibrium. These effects are much smaller than those mentioned above, and 
should be the dominant sources at low density when particle interactions are of a purely binary character. 

We present calculation of the bulk viscosity of a dilute gas from two perspectives. For each calculation we 
neglect Bose/Fermi effects and assume that relaxation can be described with a single relaxation time independent 
of the particle's momentum or species type. This is certainly unrealistic, but the single-relaxation time derivation 
allows one to precisely relate the alteration to the phase-space density, Sf{p), due to a change in the stress- 
energy tensor, STij. First, we calculate the bulk viscosity through the Kubo relation for a dilute gas. This 
is done by first finding the fiuctuation of the pressure for non-interacting particles, then multiplying by the 
relaxation time. The second calculation is based on a more dynamical picture based on the evolving phase 
space density in an isotropically expanding medium. Both calculations yield the same result. 



A. Calculating the bulk viscosity through the Kubo relation 



The Kubo relation for bulk viscosity. 



C = (/3/2) / d^x{Sf{x)Sf{0)), 
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assumes the averaging is for states with fixed energy and fixed charges. "Charges" refer to anything that is 
conserved on the time scale for which the correlation lasts. Unfortunately, it is easier to calculate thermal 
averages in the grand canonical ensemble, which allows energy and charges to fluctuate. For the GC ensemble, 
one needs to replace [T7] 
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The partial derivatives are complicated given that the pressure is most easily calculated in the GC ensemble 
as a function of the temperature f3 = 1/T and the chemical potentials for each species, = —[ig^jT. After some 
algebra, they can be expressed in terms of dfiP and da^P. Since we consider only non-interacting hadrons, with 
no Bose or Fermi statistics, the expression simplifies by using Pg — p^T, and dpP = —{P + e)/f3, 
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Assuming simple exponential decays in time, characterized by a common relaxation time, the Kubo relation can 
be calculated once one knows the equal-time fluctuations, {5A5B)^ where A and B could be any combination 
of the operators 5T, STqq or dp^. For non-interacting particles, the operators are of the form. 
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where k sums over all particles in the subvolume V. For non-interacting particles, equal-time fluctuations of 
the type {6A6B), which have sums of the type ^ are simplified by eliminating all i ^ j terms, since different 
particles are uncorrelated and are thus unlikely to be in the same subvolume. The fluctuation of T is then, 

/ d^SmSfixo = 0,x)) = 1 ^ i^, (11) 

k ^ 

with similar expressions for the other sums. One can then replace the sums with integrations over phase space 
with the phase space occupations for the species i being fe{p) = e~^^^~°''^ . In terms of the inverse temperature, 
(3, and the scaled chemical potentials, = —fi^/T, the expressions for the densities and correlators are then, 

{Pi) = I (03e"'''^"'"^^""^ (12) 

p-ae 

= {mjTKi{me/T) + 2miT^K2{mc/T)] , 

(Too) = E/ ^^E{p,me)e-^^^^^"^^^-'^' (14) 
= E {r4TK,{mt/T) + ^mlT^ K^irm/T)] 

J d^x{SToo{0)STooixo = 0,x)} = E/ ^e-^^(P''"^)-«^£(p, m,)^ (15) 

= E 2^^""' {<TK2{me/T) + mjT^Ks{me/T)} , 



j d^x{5T{0)5T{xo = 0,x)) = E j 



^'^P ^-BE(v.m,)-af P* 



(27r)3^ 9£(p,m,)2' 

= j d^x{STooiO)5Too{xo = 0,f)) - ^E™'^'^^) 



+ E I^^""'""'^^ 1 + E 7n(m,/T)2"r(-2n + 1, m./T), 

e '"^ I n=l 

7i = -l/2, 7„ = 7„_i(n - 3/2)/n 
d^p 



J d^x{6f{0)6pe{xo = 0,f)) = J ^e-^^^P'^^^-''^E{p,me) = (T,). 



(16) 



J d'x{SpeiO)dpe'{xo = 0,x)) = E/ (^^-^^(P'-^)-"^^,,,. = (p,)5,,,,, (17) 

d'x{Sf{0)6Too{xo = 0,S)) = Y.J ^"^"'^"''"'^""'y (1^) 

d-lT(dToo(0)(5Too(xo = 0,a-))/3 - E^'(P«)/3, 

y d3x(^f (0)^p,(a;o = 0,x)) = J ^e-^''^^'"'^^-"^E{p,me) = (Too,,), (19) 



(20) 



These expressions can then be inserted into an expression for the equal time commutators relevant for the bulk 
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viscosity. After lengthy manipulations using Eq. ([8]) along with the expressions for the fluctuations above, 
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where dpe = — J cPx{STqo{0)SToq{x)) is given by Eq. (15) and is related to the specific heat, dpe = —T'^dre. 
If the relaxation time At is independent of momentum and relaxation time, the Kubo relation for the bulk 
viscosity, Eq. ([s]), becomes 
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with the expectation / d^x{6T{0)6T{x)) and dfje = - J d^x{SToo{0)SToo{x)) being given in Eq.s ([T5][T6|. 



B. Alternate derivation of ^ 



The bulk viscosity can also be understood by calculating T for a gas expanding under an isotropic velocity 
gradient described by V • w for a relaxation time At without collisions, and compare it to the pressure one 
would find if equilibrium had been maintained. For the collision-less expansion, one can consider a Hubble 
expansion where the collective velocities are Vi = Xt jt. For a particle whose momentum was po at proper time 
To — yjt^ ~ 2;q, the momentum at proper time t as measured in the local rest frame, will be, in the absence of 
collisions, 



P = Po , 
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and for small time differences. 
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and for the Hubble expansion, V • w = 3/t, so that 

Ap = -p(V • w)At/3. 
One can then express the pressure and change in pressure as: 
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This last step exploited the Liousville theorem stating that dN = Vd^pf{p) / {2ttY ^ where / is the phase space 
density, stays constant. One can then use the above expression for Ap along with the fact that AV — V^(V-w)At, 
to obtain, 
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The last term in the integral can be re-expressed as an integral over energy, pdp — >■ EpdEp, and after integrating 
by parts one finds 



-(V • v)At 



-l3E(p,mi)-ai, 



(27r)3 
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with the integral being given in Eq. (16 1. To calculate the change if thermal equilibrium were maintained, one 



can apply conservation of entropy and particle number to find expressions for the changes in energy density and 
particle density during a time Ar, 



TAS = = AE ~ PAV, 
{P + e)W-v = -dte 
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(31) 



The change of the equilibrium pressure during this time is 
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Taking the difference, the effect of collisions is: 
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After applying the definition of bulk viscosity, AT = — CV • and using the expression for dP/de in Eq. 
this gives exactly the same expression for C, as in Eq. (24) of the previous subsection. 
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C. Temperature dependence of the viscosity 



According to the Kubo relations, the viscosity is a product of an equal-time fluctuation and the relaxation 
time. The relaxation time, Ar, falls roughly inversely with the density thus making both the shear and bulk 
viscosities strongly temperature dependent. The equal-time fluctuations also depend on temperature. For the 
bulk viscosity, the fluctuation, J d^x{6T{0)6T{x)) , is zero for both the non-relativistic and ultra-relativistic 
limits. This can be understood by considering the behavior of the momenta as measured in the local rest-frame, 



For an isotropic velocity gradient, 
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The corresponding change in the phase space density is proportional to 

For a change in the temperature and chemical potential, 

A/ cx -Af3Epf - Aaf, (38) 

which means that for an ultra-relativistic particle, Ep = p, the change in the phase space density is equivalent 
to a change in the temperature. The phase space occupancy is thus changed in a way consistent with changing 
the temperature, which means that the distribution remains thermal. Since a thermal distribution is unchanged 
by collisions, the collisions play no role and the coUisionless limit is indistinguishable from the equilibrated 
limit, thus there is no bulk viscosity. For a flat-space Hubble expansion one finds that the coUisionless limit is 
equivalent to having the temperature fall as 1/r. 

For the non-relativistic limit, Ep = m + /2m, and the phase space density falls in such a way that can 
be perfectly well described by changing the temperature and chemical potential. Again, this is consistent with 
local kinetic thermalization which makes the collision-less and equilibrated limits indistinguishable, and the 
bulk viscosity is zero. For the Hubble example, the temperature falls as 

However, the equivalence is broken for the case where the temperature is of the order of the mass. The 
fluctuation is then small, which leads to a small viscosity. A larger fluctuation occurs for the case of a mixture 
of relativistic and non-relativistic systems. This is most easily understood from the Hubble example above. If 
one component of the phase space density cools as while the other cools as 1/t, the two components become 
characterized by different temperatures j27) . This departure from equilibrium enables entropy production. Such 
is the case for a hadronic gas with temperatures > 100 MeV, a range for which the pion is relativistic while the 
baryons are mostly non-relativistic. 

Figure [l] displays both the fluctuation (t| defined as in Eq. 23 and the bulk viscosity, C, = /3cr^Ar, for a hadron 



gas assuming a that the relaxation time is given by the simple expression. At = l/crptot: where ptot is the total 
density of hadrons and cr is a fixed cross section of 20 mb. For this energy range cross sections are usually rather 
larger, but many of the collisions carry momentum forward which should lengthen the effective relaxation time. 
More detailed microscopic considerations are taken into account in [13j . The hadronic gas is calculated for 
two cases. First, an equilibrated gas is considered for temperatures below 170 MeV, using only the hadrons 
from the lowest lying baryon and mesons flavor decuplets and octets. Secondly, a gas is considered under the 
constraints that the net number of strange particles, baryons, effective pion number, omegas and etas are fixed 
from what the gas would have had at a temperature of 170 MeV before undergoing an isentropic expansion. In 
both cases, the fluctuations, scaled by PT, are small throughout the temperature range. At low temperatures 
the fluctuations return to zero as even the pions become non-relativistic. The two cases have very different 
densities at low temperatures. Due to the chemical non-equilibrium, the non-equilibrated state has many more 
particles due to the finite fugacities that develop [28j , which gives much shorter relaxation times and thus much 
smaller viscosities. Since the temperature range for which one would switch from hydrodynamics to microscopic 
prescriptions is in the neighborhood of 150 MeV, the behaviors shown in Fig. [TJfor lower temperatures are only 
for academic interest. For either case, the bulk viscosity is negligible in the range of interest, although it should 
be emphasized that these figures illustrate only the viscosities due to kinetic non-equilibrium. 



III. CALCULATIONS OF SHEAR VISCOSITY 



Expressions for the shear viscosity are more common in the literature and can be found in text books |29) . 
For completeness, derivations of the shear viscosity, 77, are provided here from the same two perspectives as 
the previous section. First, the shear viscosity is derived from the Kubo relation, and secondly, it is derived 
from a more kinematic perspective by considering the alteration of the phase space density in the presence of a 
velocity gradient. As in the previous section, the two expressions are identical. The latter derivation illustrates 
how one can choose a phase space density for a microscopic simulation such that it is consistent with a viscous 
hydrodynamic situation. 

Unlike the bulk case, the Kubo relation for the shear viscosity, 

r; = al ^ j d\ {ST,y{0)ST,y{r)} , (39) 



does not have the complications associated with going from the canonical to the microcanonical ensemble, 
because Tj;y is not correlated with fluctuations of the energy density or conserved charges. Thus, Eq. (39) can 
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FIG. 1: Bulk (lower left panel) and shear (lower right panel) viscosities are shown as a function of temperature for two 
cases of a hadronic gases: a chemically equilibrated gas (squares) and a second gas whose composition is that of a system 
that isentropically expanded from a temperature of 170 MeV while conserving several effective charges described in the 
text. The upper panels show the corresponding equal-time fluctuations of the stress-energy tensor elements defined in 
Eq.s Q and ( 39 1 . One reproduces the viscosities by multiplying the relaxation time and the fluctuations, then dividing 
by T. The low values of the bulk viscosity derives from the small fluctuations, which disappear in both the ultra- 
relativistic and non-relativistic limits, hence the return to zero at low temperature. For temperatures of interest for 
interfacing hydrodynamic and microscopic treatments, T ~ 150 MeV, the bulk viscosities due to kinetic non-equilibrium 
are negligible. 



be evaluated in the grand canonical ensemble. As in the previous section, the fluctuation in the dilute limit 
only comes from a particle being correlated with itself for a time Ar, which combined with the relaxation time 
leads to 

15 (27r?i)3i?(p,m,)2 



3/3Ar 



(fx{ST{0)ST{xo = 0,f)) 



with the integral being given in Eq. (16). 

One can also derive the expression for the shear in Eq. ( 40 1 from the perspective of how the phase space 
density changes in the presence of a velocity gradient. To see this one considers a particle of momentum p, 
that moves a distance {p/Ep)At between collisions. The particle then finds itself in a region whose collective 
velocity has changed by an amount, 

Ai^eoiu = ^^{pj/Ep)At. (41) 
oxj 

For small AwcoU,ij the momentum as measured in the rest frame of the new neighbors, is lessened by an amount, 
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In between collisions, the phase space density is fixed for a point in phase space following a particle's trajectory, 
i.e., Liouville's theorem. This allows one to calculate the stress-energy tensor element. 



where the Jacobian d^p/d^p' from the transformation in Eq. (42 1, p = p' + Ap, is unity. Assuming that the 
only non-zero component of the velocity gradient is dxVy, 

flA+V [ _J^LP_p'I^E(p,rn,)-a, _PxPy_ 



^\ o..^ f d^P _~8E(v>.rnA-c., PxP 

dy 

where the last step involved an integration by parts for the integral over p^- Using the definition of shear 
viscosity, one can read off an expression for 77, 



Q 2 2 
n = BAtS" ( ^ g-/?g(p,m,)-a, PxPy ,^gx 



which is the same expression derived from the Kubo relation, Eq. (40). 



IV. INTERFACING HYDRODYNAMIC AND BOLTZMANN MODULES 

It is inappropriate to model the later stages of a reaction with hydrodynamics. Once the mean free paths 
approach a similar scale to the system size, local kinetic equilibrium can be lost between different species. 
For instance, due to their lighter mass pions begin to flow outward faster than protons [30 , and due to the 
fact that they are relativistic, they also cool more slowly [27 . Although small deviations can be incorporated 
into hydrodynamics with viscosities, heat conductivities or particle number diffusions, at low temperatures, 
it becomes necessary to model with a Boltzmann prescription. The interface temperature needs to be high 
enough such that viscous hydrodynamics is warranted just above the interface temperature, while it needs to 
be low enough that hadrons are well defined objects and that the dynamics are not greatly affect by non-binary 
interactions. Albeit narrow, there does seem to exist such a window with 140 ^ 165 MeV. 

In order to initialize the Boltzmann description, one needs to choose a phase space density, /(p,a;). Since 
/(p, x) depends on the three-dimensional momentum p at any space-time point x, it can be altered in innumer- 
able ways to incorporate viscous corrections, so long as one maintains a consistency between the stress-energy 
tensor elements of the input hydrodynamic model and that of the microscopic model, 

species t 

where £ refers to particular species (spin and flavor degeneracy factors are suppressed). The functional form 
for the modification of / might or might not be based on a physical picture. A common ansatz suggested by 
Grad 3lj (motivated by the expansion of the phase space density in a viscous fluid) for correcting for shear is 

/(P, X) = /eq(p, X) [1 + Cip)p,Pj7Tlf] , (48) 

where /eqis the equilibrated phase space density determined by the energy and charge densities. As discussed 
in [5], the momentum dependence of C{p) can be picked to reflect the energy dependence for which particles 
lose memory regarding their original momentum. Additionally, it might depend on species type, e.g., since 
protons have larger cross sections than pions, they have relatively smaller viscous corrections to their phase 
space density. 

The advantage of Eq. [48] is that it is related in a straight-forward manner to the viscous part of the stress- 
energy tensor (whether calculated in Navier-Stokes or the more general, dynamical, Israel-Stewart), as it is 
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essentially an estimate of the Grad correction term. A disadvantage of this form in Eq. (48) is that for high 
momentum and for certain directions of p one finds negative phase space densities. The unphysical negative 
phase space density is not surprising: The typical departure of thermalization varies with the momentum of the 
particle (typically, for p ^ T a ^ so At ^ p^). Thus, for large momenta, the first term of the gradient 
expansion will be insufficient. 

For the purpose of this work, we shall consider an alternative prescription, which, as we will see, is less affected 
by the high-momentum pathological behaviour when implemented in a Monte-Carlo code. In this prescription, 

f{p)=exp{-E{p')/r -a'}. (49) 

The primes on the chemical potential and temperature note that the temperatures and chemical potentials can 
be different than the quantities one would choose to match the energy and charge densities in a non-viscous 
theory. The momenta p and p' will be related to one another through the relation, 

Pl^Pi+ KjP'y (50) 



Note that for Xijp'j ^ 1, ie for small momenta, the first Taylor coefficient of this formula should lead to Eq. 48 
Only for bulk effects, i.e., J^ii^^a ~ P) ¥^ 0, will the temperature T' and chemical potential a' differ from lEe 

equilibrated values. An advantage of this form is that it can be applied to Monte Carlo generation of momenta 

in a straight-forward manner by generating p' according to a thermal distribution, then performing the linear 

transformation to generate the momenta p. 

To motivate the linear behavior in Eq. (50), one can consider a particle with momentum p' at time tq. After a 

small time At which can be thought of as the collision time or the relaxation time, the particle will have moved 

a distance Ar = {p' / E{p'))At. At that time, it will have moved to a region with velocity Avi — {dvi/drj)Arj. 

The momentum as measured in the new frame will be: 

p,^p[~ E{p')Av., = p', i?(p')|^^Ar =p',- ATd,v,p'^, (51) 

and 

X,j = -ATdjV,. (52) 

For irrotational flow, the Navier Stokes equation assumes that the velocity gradient is proportional to the viscous 
correction to the stress-energy tensor. In particular, we will assume 

A,, = ^(«)7rl;) + A(^)^W<5,„ (53) 

and will find the coefficients A'^^'^ and A'^''^ that satisfy the consistency requirements of Eq. (47) in the limit 
that A, or equivalently tt^''^ and TrC*), are small. 

With this transformation, one can apply the Liouville theorem, /(p) = /(p')) 

where it has been assumed that the volume is fixed. The stress-energy tensor changes due to A and due to 
AT = T' -T and Aa = a' - a. 



^ J (27r)3 \ 3E{p') 3E{p'r 

+ J (27r)3" 3£;(p') I 

Here, A^**^ is the traceless part of A and A^**^ describes the remainder, i.e., 



A.,-A|f + AW<S,,. (56) 
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Since the modifications must be made at fixed energy and particle densities, one also has the constraints 



(2^ 



-i5(p')/T 



-Aa( 



Eip')AT 
2^2 



r2 £;(p') 



(57) 
(58) 



Using Eq. (53), one can sub stitu te for A*-*-* and A*-*"^ in the expressions for tt^'*^ and Tre-) in Eq. (^55^ and after 
using the constraints in Eq. (57) solve for A'''^^ and A''''\ After some lengthy algebra, 



AT 

Aai 



At 

At 

3C ' 

D ' 
_AWPe, 

Dp, 



(59) 



where D is given in Eq. (|9|, the ratios ry/r and C/t are given by Eq.s (461 and Eq. (24). The results of Eq.s 
([59|) are consistent with expectations from the Navier Stokes equation. If the expressions for A are replaced 



with the product of the relaxation time and velocity gradient from Eq. (52), then Eq. (53) can be equated to 
the Navier-Stokes equation once Eq.s ( 59 ) are applied. 



The method to consistently generate a phase space density for a thermal distribution at fixed energy densities, 
particle densities and at a given departure from equilibrium of the stress-energy tensor, Tr^-, can be summarized 
as follows. First, one finds the temperature and chemical potentials corresponding to equilibrated energy density 
and particle densities. If tt*^''-' is non-zero, one must then alter the temperature and chemical potential according 



to Eq. (|59|) . After generating a thermal distribution, each momentum is scaled by the matrix Sij + Xij , where 



Xij is given by Eq.s (59) and (53) 



The work above demonstrates that generating an equilibrated phase-space density, followed by a rescaling of 
the momenta using matrices built proportional to tt^^^ or 7r'''\ one will generate stress-energy tensors consistent 
with the TT^**^ and tt'^''^ used to scale the momenta. Further, the constants of proportionality relating A and 
TT can be found by integrating moments of the the thermal distributions, given in Eq.s (53). However, these 



relations were derived by keeping only the linear terms in A when calculating the stress-energy tensor in Eq. 
(55). The consistency from the above prescription should break down for large enough A, or equivalently large 
enough tt'-"-' or tt^''^, as the relation between A and tt will no longer be purely linear. To see the range to which 
the linear behavior extends, we consider an equilibrated gas of hadrons at a temperature of 160 MeV, with the 
phase space density modified according to the following form for A, 



X — - X - X — 4(s)„('"P"*) 

^xx ^yy ^mag ^ ; 

with all other Xij = 0, and oi references one of the deviations to the stress energy tensor. 



as 
b 



{Txx ^ 2^yj/)/2, 0,2 



1 



12 



T 



(60) 



(61) 



T 

' T.' 



(24 



T 



P)/3. 







For the procedure outlined above to be successful the stress energy tensor generated with the specified A^ 
should result in ai — a^'"'^^*'', where a^™''^'"' refers to input value from the hydrodynamic simulation, whereas 
ai refers to result from integrating the phase space distribution. For the test, it is assumed that a^^^j^"* 

and 6('"P"*) = 0, so that a^^i, h and Ae should all be zero. For small a^'"^"'^ the procedure should be exact, 
but for larger deviations non-linear contributions will lead to the output values varying from the input values. 
This is illustrated in Fig. [i] which shows a^, b and the change in the energy density Ae as a function a^"^^"*^ 
Given that systematic uncertainties, both theoretical and experimental, at RHIC are rarely below the 5% level, 
any deviation of the order of one percent or less should be tolerable. As can be seen in Fig. |2j the procedure 
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0.1 0.2 0.3 0.4 0.5 

a^/P (input) 



FIG. 2: (color online) For the procedure to reproduce a phase space density consistent with the stress energy density 
and particle densities, momenta were generated according to a thermalized distribution, then scaled linearly according 
to a matrix Xij oc TTij, and with the coefficient of proportionality chosen so that the procedure is exact for small input 
values, tt'"^''"'''. Here, the resulting values of the stress-energy tensor, as defined in Eq. (61 1, are shown as a function 
of the input values aj^'"^"*^ which could come from a hydrodynamic simulation. For the test, the other input values are 
set to zero, i.e., ^('"P"*) — Sj"''"'' = 0. The resulting value of ai (circles) closely follows the input value (line) over the 
same range. The other elements, b (upward triangles), a2 (green squares) and Ae (downward triangles) stay near zero 
as desired for small aj^'"^"*', with relative deviations not staying below one percent for a5^'"^"''/P < 1/3. 



described here appears sufRcient for values of a^™''"*''/P ^ 1/3. Although not shown here, similarly acceptable 

violations of the particle densities occur for large a^'"''"*''/P ^ 1/3. For larger deviations, one would have to 
either accept the discontinuities of the densities and of the stress-energy tensor, including the energy density, 
or apply a more sophisticated, and likely more difficult to implement, prescription. 

Along the same lines as the tests illustrated in Fig. [2] the validity of the linear considerations for correcting 
for bulk effects was also analyzed. Due to the small bulk viscosity for a gas, the input values of 6('"p^') /P 
are small, but nonetheless can require fairly large modifications of the phase space density since the linear 
coefiicient A''') w l/C- bmce are not small, non-linear contributions tend to have similar corrections to 
the stress energy tensor as was the case for the shear viscosity. However, since the linear corrections are so 
small, the non-linear corrections quickly overwhelm the linear corrections and invalidate the prescription. Thus, 
while the linearity assumed in the procedure worked well for generating a consistent stress-energy tensor for 
the shear case even for fairly large values of a'"™^^^^ /P, the linear procedure failed for the equivalent bulk case. 
For instance, if one were to set the viscosities according to the Navier-Stokes equation in a one-dimensional 
expansion, the magnitude of the components of A would be of the order of the relaxation time multiplied by the 
velocity gradient, for both the shear and bulk corrections. The linear procedure for correcting for both the shear 
and bulk terms would both lead to unwanted non-linear deviations of the energy and particle densities, and to 
the stress energy tensor. These deviations would be of the same order for both the shear and bulk corrections. 
However, for typical values of A seen in heavy ion collisions, the unwanted deviations to the stress-energy tensor 
tend to be much smaller than the input value a''"^"*'', while they tend to be much larger than 5('"P"*^. Thus, 
there is little point to correcting for the bulk effects since the induced errors tend to be larger than the error 
associated with simply setting ^('"P"') = 0. 

To generate particles according to the phase-space density and the evolution of the breakup hyper surface, 
one can apply the relation [5], 

dN= (2^^j'i^(p) (P-")Q(P-")' ^"^^''^'''^XpAY^AZs. (62) 

The vectors AX, AY and AZ describe the widths of the "volume" element. For a surface element of the 
breakup hypersurface, defined by a quadraleteral, one can consider the corresponding surface element along 
the hypersurface at a later time. The three vectors describe the separation between the centers of each of the 
opposite faces. To understand how this provides a volume element, one can consider the hyper surface element 
defined simultaneously across the element, i.e., AX^ — AY^ — 0. The product e"^''^ AXjAYs then becomes 
an antisymmetric matrix with the "iO" elements defining a vector perpendicular to the surface element whose 
magnitude is the area AA. If AZ^ = 0, the boundary of the next surface element for the breakup surface is 
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defined simultaneously with the previous one, and the product t'^^'^^ ^Xf^^y^^Z^ is purely time-like with a 
magnitu de equal to a volume Ay bounded by the concurrent hyper-surface elements. The final convolution 



in Eq. (62) thus gives dN — fAV(Pp/{2TT)^, as desired. If the hyper-surface element is defined at the same 



position, but different time, i.e., AZ — (At, 0,0,0), the procedure yields dN — /AA -pd^p/ep, which is what 



one expects for emission from a static surface. The step function Q{p ■ fl) in Eq. 



62| comes into play only 



for space-like il" , and eliminates those particles traveling inward into the surface [34|, I35j [37] . For collisions at 
RHIC this amounts to only a few percent of all the particles. The step function must still be treated carefully, 
as the sign can depend on how AX and AY are defined. For space-like Afl, one should ensure that one is 
choosing those particles leaving, rather than entering, the breakup surface. 

The algorithm for creating a sampling of particles from the element Af2 using the description above is: 

1. Calculate the "volume" vector Afl" as described above. 

2. Assuming a real volume of size V = \Ail\, calculate the number of particles N one would create in a 
thermal system with that volume, temperature and chemical potential. 

3. For each particle in N, create a momentum consistent with a static thermal distribution. 

4. Scale the momentum according to Xij determined from from the deviations of the stress-energy tensor as 



described in Eq. ( 53 1 . 



5. Boost the particle according to the velocity of the fluid element. 

6. If p • O < 0, throw out the particle. One must be careful to ensure that p-^l > refers to particles leaving 
the element. 

7. Keep or reject the particle according to the probability p ■ fl/EpV. 

It would be straight-forward to alter the procedure above to incorporate non-uniform relaxation times, i.e. A 
could be a function of the magnitude of p as measured in the fluid frame or could be a function of the species. 
The principal motivations for altering the form would be to incorporate longer relaxation times for high energy 
particles, given that they may require multiple collisions to re-thermalize, or to give longer relaxation times 
to species with longer mean free paths. It has been shown that the difference between saturating V2 values 
for different species, often associated with quark-number scaling in coalescence, might partially derive from 
different relaxation times j5] . If two species are emitted from the same fireball, but characterized by different 
relaxation times, the species with the shorter relaxation time (e.g. protons) will have a more isotropic local 
momentum distribution. Since the anisotropy of the local phase space density lowers V2, those particles with 
smaller relaxation times will then have higher values of V2 ■ 



Comparing this procedure to the Grad form described by Eq. ( 48 1 , reveals no difference if the departure from 
equilibrium is small, assuming an appropriate choice is made iorC(p) in Eq. (48 1. To see this, one can expand 
the argument in the exponential for the phase space density for small A, 

/(p') = e-^(P')/^ « e-^(P)/^ {l + ^A,;,,} (63) 



where the definition for A^*' in Eq. (53 1 was used to replace A. By comparison with Eq. (48) one can see that 



the two expressions are identical for small tt if 



Cip) - (64) 



EpT 



For the assumption of uniform relaxation time applied thus far for the scaling procedure, A^''' was indpedendent 
of the momentum. However, one could in principle choose an arbitrary function form, i.e., instead of the 



assumption that the scaling function A is independent of momentum as in Eq. ( 53 1 , one could instead have 
chosen, 

A,,(p) = F(p)4^V,„ (65) 



where F{p) could be any arbitrary function of the magnitude of the momentum (recent developments show 
F{p) could indeed be non-trivial and species-dependent [35]) One would then retrace the procedure above and 
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derive the value for Ap\ In the smah tt^^^ limit, the Grad form and the scaling procedure would then equate 
once one had chosen 

C{p) = ^ y ■ (66) 

Even if one were to use the Grad form rather than the scaling procedure presented here, this expression offers 
insight as to how one might choose the form for C{p). If one thought the relaxation time were independent 
of momentum, or particle species, one would choose C{p) consistent with F{p) being constant. Otherwise, it 
would make most sense to choose F(p) proportional to the relaxation time. For higher momentum hadrons, 
one might expect longer relaxation times, not so much because of reductions in the overall cross sections, but 
because it can take several collisions for a high energy particle to completely lose memory of its momentum. 
The extreme case would be a high-energy jet, for which the relaxation time might be so large that it cannot 
re-thermalize on the time scale of the expansion. In this limit both the Grad form and the scaling procedure 
lose their validity. 



The Grad form, Eq. ( 48 ) , has drawbacks in that it yields negative phase space densities for large momentum. 
The momentum range for which the phase space densities become unphysical depends on the size of the viscous 
corrections. For correction of the order oi/P ~ 1/3, negative phase space densities ensue for momenta only a 
few times the thermal momentum. This problem occurs even for the case where F{p) is chosen to be a constant. 
If one chooses a functional form F{p) in Eq. ( |66[ ) that grows with momentum, the problems of negative phase 
space densities become all the more acute. 

The scaling algorithm outlined above also has some drawbacks in that one loses the equivalence between the 
non-equilibrium components of the generated stress-energy tensor and the non-equilibrium components used 



to determine the coefficients in Eq. (53) is lost for large departures from equilibrium. One advantage of the 



procedure described above is that it accommodates Monte Carlo generation of particles in a fairly straight- 



forward manner. The procedure is also better behaved at large momentum. From Eq. (65), one can see that 
if all particles are given a fixed lifefime, i.e. F{p) is a constant, that the scaling matrix A is independent of 
momentum, and as long as A is not large, the linear approximation is reasonable for all momenta. However, 
if F(p) increases with momentum, there will again be some momentum range for which the approximation 
becomes unwarranted and perhaps even unphysical. 

Both procedures have difficulty adjusting for bulk corrections to the stress-energy tensor, i.e. ^ 0, 

as one must also adjust the temperature and chemical potentials to maintain fixed energy and particle number 
densities. Once these densities are conserved, both forms require large deviations of the phase space density 
to produce small h/ P. Although it would be straight-forward to do this consistently, linear approximations 
(assuming that the modification factors are proportional to tt^''-'), are bound to run into difficulties such as 
negative phase space densities for specific momenta in the Grad form, or for non-linear corrections overwhelming 
the linear corrections for the scaling procedure. Given that the bulk corrections to the stress energy tensor tend 
to be less than one percent of the pressure in the region of hydrodynamic/microscopic interface, the best course 
of action may be to ignore the bulk correction altogether. This recommendation would change if the bulk 
correction were not small for the hadron case. Such might be the case if the hadronic mean fields |17) . e.g. the 
chiral condensate, are away from equilibrium. However, that would suggest adding a dynamic mean field to 
both the hydrodynamic and microscopic prescriptions [26', rather than attempting to incorporate the viscous 
corrections into the momentum distributions, which as shown in Sec. |ll| can accommodate only small viscous 
corrections. 



V. SUMMARY 



By calculating the bulk and shear from both a dynamic perspective and from Kubo relations, it is clear how 
one should alter the phase space distribution for a system where the relaxation time is independent of species 
or of momentum. The calculations of the bulk viscosity were especially useful, as they illustrate why the bulk 
viscosity is so small for most applications to heavy ion collisions. In particular, one can see that C — ?> for 
both ultra-relativistic and non-relativistic system. A mixture of the massive (T << m) and massless (m << T) 
components leads to the largest values of C/P-: but even in those cases bulk viscosity from kinetic non-equilibrium 
is negligible for designing an interface between hydrodynamic and Boltzmann models in heavy ion collisions. 

The viscous calculations alluded to above inspired the algorithm described in Sec. |IV| for interfacing hydrody- 
namic and microscopic descriptions. The method was based on generating particles with a thermal distribution, 
then scaling their momenta linearly with a matrix composed of a unit matrix plus a piece proportional to the 
non-equilibrium stress-energy tensor taken from the hydrodynamic calculation. The constant of proportionality 
was then expressed in terms of integrals of the thermal distributions. The method was shown to be better than 



15 



one percent accurate as long as the non-equilibrium shear components of the stress-energy tensor remain ^ one 
third of the equilibrium pressure. The method did not appear useful for accurately generating non-equilibrium 
deviations of the bulk pressure, although such deviations are so small they can be safely neglected. The method 
has an advantage over previously considered techniques to generate corrected phase space densities, in that the 
phase space densities remain positive for all momenta. The proposed method also lends itself well to generating 
Monte Carlo ensembles of discrete momenta. 
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